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Abstract — This paper starts by an investigation of nonlinear 
transmission in space-division multiplexed (SDM) systems using 
multimode fibers exhibiting a rapidly varying birefringence. A 
primary objective is to generalize the Manakov equations, well 
known in the case of single-mode fibers. We first investigate a 
reference case where linear coupling among the spatial modes of 
the fiber is weak and after averaging over birefringence fiuctu- 
ations, we obtain new Manakov equations for multimode fibers. 
Such an averaging reduces the number of intermodal nonlinear 
terms drastically since all four-wave-mixing terms average out. 
Cross-phase modulation terms still affect multimode transmission 
but their effectiveness is reduced. We then verify the accuracy of 
our new Manakov equations by transmitting multiple 114-Gb/s 
bit streams in the PDM-QPSK format over different modes of a 
multimode fiber and comparing the numerical results with those 
obtained by solving the full stochastic equation. The agreement 
is excellent in all cases studied. A great benefit of the new 
equations is to reduce the computation time by a factor of 10 or 
more. Another important feature observed is that birefringence 
fluctuations improve system performance by reducing the impact 
of fiber nonlinearities. The extent of improvement depends on 
the fiber design and how many spatial modes are used for SDM 
transmission. We investigated both step-index and graded-index 
multimode fibers and considered up to three spatial modes in each 
case. In this paper, we also consider multicore fibers and show 
that the Manakov equations remain valid as long as the length 
scale of birefringence fluctuations is about 100 times shorter 
than that of linear coupling. Finally multimode fibers with strong 
random coupling among all spatial modes are considered. Linear 
coupling is modeled using the random matrix theory approach. 
We derive new Manakov equations for multimode fibers in that 
regime and show that such fibers can perform better than single- 
modes fiber for large number of propagating spatial modes. 

Index Terms — Birefringence, Manakov equation, Multimode 
Fiber, Multicore Fiber, Space-division multiplexing. Fiber Non- 
linearity. 

I. Introduction 

Recent work on the optical fiber capacity Hmit [l]-[3] has 
indicated that the continued use of single-mode fibers for 
modern telecommunication systems may not be able to support 
the growing data-traffic demand in the near future [4], [5]. 
Although more single-mode fibers can be deployed as a short- 
term solution, it is important to look for alternative solutions 
for the next generation of optical transmission systems. The 
impressive bit-rate increase of the last two decades was 
made possible by exploiting diverse properties of the elec- 
tromagnetic field through wavelength-division multiplexing 
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(WDM) in combination with phase modulation formats and 
polarization-division multiplexing (PDM). Recently, multi- 
mode fibers have become the focus of attention [6], [7] as they 
permit simultaneous propagation of multiple spatial modes. 
Multiplexing signals using different spatial modes of an optical 
fiber is referred as space-division multiplexing (SDM). Recent 
work on multimode fiber has focused on demonstrating the 
feasibility of the SDM technique and how digital signal 
processing (DSP) can help to compensate linear distortions of 
the signal [8], [9]. However, nonlinear propagation effects are 
the principal limitation of optical fiber communication systems 
and their understanding in the case of multimode fibers is still 
very limited compared to the case of single-mode fibers. 

A set of nonlinear propagation equations has been derived 
for multimode fibers in [10] by expanding the optical field 
in terms of vectorial fiber modes. It is difficult to use these 
equations in their most general form to study birefringence 
effects. In the case of silica fibers used for telecommunication 
systems, we can simplify these equations by making use of 
the so-called weakly-guiding approximation. In this paper 
we derive in section II a new set of multimode nonlinear 
equations, which represents the two polarization components 
of each spatial mode through a Jones vector and includes 
linear coupling between spatial modes. We then use them to 
study the effects of fiber's random modal birefringence and 
linear mode coupling using random matrices in two important 
cases of practical interest that we refer to as the weak and 
strong coupling regimes. In the weak-coupling regime, linear 
coupling among distinct spatial modes is weak compared to the 
birefringence-induced coupling between the two polarization 
components of the same spatial mode. In the strong-coupling 
regime, both types of coupling occur in a random fashion with 
the same order of magnitude. In practice, certain spatial modes 
of multimode fibers may be weakly coupled whereas others, 
may be more strongly coupled. For instance, it was observed in 
[8] that the LPl la spatial mode is strongly coupled with LPl lb 
but weakly coupled with LPOl. The strong and weak spatial 
mode coupling regimes studied in this work represent two 
extreme cases for practical systems. Multimode fibers always 
have, in general, some degree of coupling between any pair of 
spatial modes. Important factors impacting how much linear 
coupling there is between spatial modes are: fiber design and 
fabrication, cabling, bending and twisting. 

It is well known in the case of single-mode fibers that 
rapidly varying birefringence results in an averaging effect 
that reduces the strength of nonlinearities: mathematically, 
this averaging is described by the Manakov equations [11]- 
[13]. Its use also reduces the computation time considerably 
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when the nonUnear propagation is performed numerically. In 
general, the standard Manakov equations cannot be used in 
the case of multimode fibers since they only describe how 
birefringence impacts intramodal nonlinear effects and do not 
include intermodal nonlinear effects. In section III, we derive 
a new set of Manakov equations for nonlinear propagation 
in multimode fibers in the regime of weak coupling between 
spatial modes, following the ideas in [11] and averaging over 
all possible polarization states. In section IV, the impact of 
rapidly varying birefringence on a specific SDM transmission 
system is studied through numerical simulations. We show that 
the new set of Manakov equations successfully predict the 
effects of random modal birefringence, while also reducing 
the computational complexity of the problem. In section V, 
we apply the Manakov equations to the case of a multicore 
fiber whose cores are linearly coupled with a constant coupling 
coefficient. We show that Manakov equations remain valid 
in the weak coupling regime as long as the length scale of 
the birefringence fluctuations is 100 times shorter than that of 
linear coupling. 

Finally, we consider in section VI a case with large random 
coupling between spatial modes, referred to as the strong- 
coupling regime. Although current transmission experiments 
over multimode fibers are not fully in the strong-coupling 
regime [6], this case has attracted theoretical attention to 
obtain an understanding of stochastic effects such as ran- 
dom differential group delays among various modes [14] 
or random mode-dependent loss within the multimode fiber 
[15]. Recently Manakov equation have been generalized for 
multimode fibers in the high coupling regime [16]. Following 
a different approach based on random matrix theory, we obtain 
similar results within the weakly guiding approximation. A 
comparison of the weak- and strong-coupling regimes is also 
presented in this paper. 

II. Nonlinear propagation in multimode fibers 

We simplify the analysis in [10] by making the weakly 
guiding approximation, known to be quite accurate for silica 
fibers with a relatively low difference in the refractive indices 
between core and cladding. In this approximation, the modal 
fields are linearly polarized in the transverse plane, as the field 
components along the direction of propagation (z axis) are as- 
sumed to be negligible. Moreover, the spatial modes associated 
with the two polarization components are assumed to have the 
same spatial distribution. With these simplifications, we write 
the total electric field as a sum over the M distinct spatial 
modes (we use the same definition of the fiber modes as in 
[8]) of the fiber in the spectral domain: 

^ M ^ 

E(x, y, z,uj) = Y, e^^-(-)- A^(z, c^)F^(x, ^)/ V^, (D 

m 

where A^(z,cc;) = [A^x(^,^) A^y(z,cj)]^ is the Fourier 
transform of the time-domain slowly varying field envelope 
of the mth mode represented in the Jones vector notation. 
The superscript ^ denotes the transpose operation and the 
tilde denotes a frequency-domain variable. It includes the 
slowly varying amplitudes of both polarization components 



for the spatial mode m with the spatial distribution F^(x,?/) 
and the propagation constant (5m{oo), expressed in the form 
of a diagonal matrix to account for fiber birefringence, i.e., 
(3m = diag[/3mx Pmy\\ 6*^"^^ corrcsponds to the usual defi- 
nition of an exponential matrix. Note that ¥m{x^y) is real in 
the weakly guiding approximation, and the phase of the optical 
envelope is contained in the complex field A^. Finally, the 
normalization constant in (1) represents the power carried 
by the mth mode. It can be expressed as = ^ ^eff cim, 
where Im = fim/nei^ JJ F'^{x^y)dx dy, cq is the vacuum 
permittivity, neff is the fiber effective index of the fundamental 
mode and flm the fiber effective index of the mode m. 

Each frequency component of the optical field satisfies the 
frequency-domain Maxwell equations, 

V^E{uj)^n\x,y,z)^E{uj) = -c^VoP^^M, (2) 

where n{x^y^z) is the refractive index distribution within the 
fiber whose z-dependence accounts for small perturbations 
along the propagation distance, i.e., n^{x^y^z) = nKx^y) + 
An'^{x,y^z) [17]. P^^(a;) is the Fourier transform of the 
third-order nonlinear response [18] 

pNL(^) = jX^^^ ([E^EJE* + 2[E^E]E) . (3) 

Here * denotes the complex conjugate and the superscript 
^ denotes the Hermitian conjugate, x*^^^ is the instantaneous 
nonlinear response of silica. 

We assume that the spatial distribution of a mode is not 
perturbed by nonlinearities or by small variations of the 
refractive index along the propagation distance z. In that case, 
Fjn{x^y) satisfies the eigenvalue equation 

\/^Fm{x,y) + (^^nl{x,y) - Fm{x,y) = 0, (4) 

where the subscript i = x, y refers to the two orthogonal 
polarizations of the mth spatial mode. In fibers with per- 
fect cylindrically symmetric refractive index profile no{x^y), 
z^mx = Pmy In practicc, because of imperfections, fibers 
exhibit some polarization birefringence and /3^x 7^ Pmy 
In the remainder of this section we consider a sufficiently 
short fiber section such that the birefringence remains constant 
along the fiber length. The case of accumulation of randomly 
oriented birefringent fibers segments is treated in subsequent 
sections. Note that the various solutions of (4) are chosen to 
be orthogonal to each other and are normalized such that 

J J Fm{x,y)Fp{x,y)dxdy = lmSmp^ (5) 

Smp being the Kronecker function. 

To isolate the evolution of a specific spatial mode, say the 
pth mode, we use their orthogonal nature. We substitute (1) 
in (2), multiply both sides with Fp, and integrate over the 
transverse x-y plane. The resulting equation is then converted 
to the time domain by following a standard procedure and ex- 
panding (5p{(jj) in a Taylor series around the carrier frequency 
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ujQ [19]. The final result can be written in the form 

dAp I .dAp l32pd^Ap 

^= ^{^5,p-^)Ap - - 



Imn 



QmpAf^ 



(6) 



where Ap{z,t) is the time-domain slowly varying envelope 
of the pth mode expressed in a reference moving frame at 
the group velocity Vg^. Pr is an arbitrary reference propa- 
gation constant for all the spatial modes and /3op = /3p(co'o), 
l^ip = d(3p/duj\ujo, hp = d'^(3p/duj'^\oj^ are respectively 
the propagation constant, inverse group velocity and group- 
velocity dispersion (GVD) of the pth spatial mode. We assume 
here that the two polarization components of a spatial mode 
may have different group velocities but have the same GVD. 

The nonlinear parameter in (6) is defined in a similar fashion 
as in single-mode fibers, by 7 = 3/coX*^^V(4eocng^Aefr), 
where Agff is the effective area of the fundamental mode. 
The linear and nonlinear coupling among spatial modes is 
governed, respectively, by 

qmp{z) = — — ,^ II An^{x,y,z)F^Fpdxdy,(7) 
i^p ) J J 



flmnp - 



2neff(Iri 
Aeff 



{Idmlnlpy^^ 



FiFmFnFp dx dy. 



(8) 



Note that linear coupling terms appear as a result of the 
perturbation of the refractive index An^(x,?/,z). In the case 
of an ideal fiber, An^(x, y^ z) = 0, no linear coupling occurs 
between the spatial modes. 

Equation (6) governs propagation of arbitrarily polarized 
light in the spatial mode p within the Jones-matrix formalism. 
It includes all third-order nonlinear effects (both intramodal 
and intermodal kinds) as well as dispersive and birefringence 
effects. Fiber losses can also be included by adding the term 
— {ap/2)Ap on its right side, where the loss parameter ap can 
be different for different modes. 

III. Manakov Equation in Weak-Coupling Regime 

Equation (6) assumes constant birefringence along the fiber 
length. In practice, owing to fiber imperfections such as a 
nonuniform core whose shape varies along the fiber, bire- 
fringence varies randomly and rapidly on a length scale that 
is short compared to the effective lengths associated with 
the GVD and various nonlinearities. This feature can be 
implemented numerically by rotating the principal axes of the 
fiber periodically after a distance shorter than the length scale 
of birefringence fluctuations. However, a numerical solution 
of Eq. 6 then requires a relatively small step size and is quite 
time-consuming. 

In the case of single-mode fibers, the problem has been 
solved by adopting the well-known Manakov equation [11] 
[12]. The idea is that, as a rapidly varying birefringence 
changes continually the state of polarization (SOP) of propa- 
gating light in a random fashion, one can average the prop- 
agation equation itself over all polarization states. We can 



follow the same procedure for multimode fibers by assuming 
that the SOP of each spatial mode evolves, randomly and 
independently, i.e., we assume that the SOP of each mode 
evolves in an independent fashion. This can be justified by 
noting that, as the spatial distributions of the modes are 
different, the influence of local stress and fiber imperfections 
may also be different from one spatial mode to another. 

A rapidly varying fiber birefringence can be tracked using 
the transformation 



Ap{z) = Yip{z)Ap{z) 



(9) 



where Yip{z), a unitary matrix belonging to the group U(2) is 
a Jones matrix of the form 



R^Rp — I2 



(10) 



(11) 



where I2 is the 2 x 2 identity matrix and Vijp are random 
variables. An important issue is how to construct Rp to ensure 
that the SOP probability for the pth mode is uniformly dis- 
tributed over all possible SOPs. We make use of the following 
procedure applicable to the general case of nxn matrices 
[26]: pick a matrix ]V[ whose all elements are normally 
distributed (with zero mean) and apply a QR decomposition, 
]V[ = RT, such that T is an upper triangular matrix with 
positive diagonal entries. Then, R is a Haar matrix whose 
elements are uniformly distributed over U(n). 

The use of (9) in (6), leads to the following equation for 

Ar)\ 



dA^ 



dAp f32p d^Ap 



dz -^^p-p -^^p dt 2 dt^ 

-^"f^^^ flmnp- {[AJ Kjll^Am] R^R*A* 



- 2 [Ap Rp R^A^] R^R^A^ 



(12) 



where 



(9R 

S^op = Rp (/3op - Rp - zR^ ^ , (13) 
5/3i^ = R^(/3i^-— )R^, (14) 

Qmp — ^mpRp Rm* (1^) 

Equation (12) is stochastic because Rm(^) is a random 
matrix that changes along the fiber length on a length scale 
associated with birefringence fluctuations. As a result, the 
birefringence parameters appearing there vary randomly. In 
addition, the intramodal and intermodal nonlinear couplings, 
as well as the linear coupling, also become random. 

In the following analysis, we average the propagation equa- 
tion (12) over all possible realizations of the matrix R^. 
Averaging over birefringence amounts to assuming that other 
phenomena that produce z-dependent variations, in particular 
the linear coupling between spatial modes (qmp{z) in (15)), 
occur over a much longer length scale in comparison to the 
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length scale of polarization fluctuations within each spatial 
mode in the weak-coupling regime considered here. 

Let us focus first on the nonlinear terms in (12) by writing 
them as 



Imn 

- 2 [ApRpR^A^] R^R^A^ 



(16) 



The matrix multiplications appearing there can be expressed 
as 



R/^R^ 



ml Im- 



RpR^ — 



aim h 

m 

/)* r 

.^ml Hm\ 



where we have defined the following quantities: 



aim 


* 

=^111 


^llm " 


* 

r ^211 ^21m 


(18) 


hm 


* 

=^111 


^12m " 


* 

r ^211 ^22m 


(19) 


^Im 


* 

— ^211 


^21m ~ 


r ^221 ^22m 


(20) 


aim 


= ^111 


^llm " 


1" ^211 ^21m 


(21) 


hm 


= ^111 


^12m " 


^ ^211 ^22m 


(22) 


^Im 


= ^211 


^21m ~ 


1" ^221 ^22m- 


(23) 



Using (17) in (16), the nonlinear terms for the x-polarized 
pth mode become: 



A * A A I ^i* A * A A 

-^pn'^lx'^mx'^nx ~r ^ml^pn-^ly^my^-^ny. 



Imn 
\^lm(^: 

A-h n A*A A ^ r n A*A A 

'^lmpn l:K my nx ' ^Im^pn ly my ny. 

I 7 * A * A A _L ?i* ?i* A * A A 
"T^Zm^pn'^Zx'^mx'^ny ^ml^pn^ly^myi-^ny 

A_ h 7) A* A A \p ^* A* A A 
' ^Im^pn ly. my ny ^ ^Im^pn ly my ny 



n n* K K A* A A A* 

^Im^pn-^lx-^my^ny ^ ^ml^pn-^ly^my^nx 

-\-h n* A A A* ^ r /7* A A A* 

' ^Im pn ly my ny ' ^Im pn ly my ny 

-\- n h* A A A*^/)*/)*AA A* 
'^^ImKn^lx^my^ny + ^Im^n^ly^my^ny ) 

(24) 

We are interested in finding the ensemble average, or the 
expectation value, that we denote by ( ) of the various 
coefficients appearing in (24) under random unitary transfor- 
mations. When the subscripts m,p, n correspond to different 
spatial modes (I ^ m p n), the coefficient aimapn 
is a sum of products of functions rn/, rn^, run, and 
^2iz, ^2im, ^2in, ^2ip, all of which are independent random 
variables with zero mean. As a result, {aimapn) = 0. With 
the same kind of reasoning, we find that all terms in (24) 
average to zero, except those containing the combinations 

^ (25) 



(bppbpp) — 1 {appCLpp)-- ^ 



^mm^pp ~ ^mm^PP ~ appttpp — Cppttpp — l 



i^pm^pm) {^pm^pm} ' 



(26) 
(27) 



i^mp^m) - {KmKm)-^' ^^^^ 



These equations are obtained using the properties of Haar 
matrices, which are given by Eq. (37) to (39) in the case of 
2M X 2M matrices and assuming that R^ are independent 
random unitary matrices. Finally, by gathering the nonzero 
terms together in ATpx and JVpy , we obtain the following 
averaged propagation equation: 



dAp 

dz 



-{S^Op)Ap 



— ^1 (^fpPPP g I Ap I 



m^p 



I32p d^Ap 



2 df2 
-lA P 



(17) with 



{Sf3op) 
{Sf3ip 



) = o( 



1 

2 
1 



Ppx 



' Ppy) 

Ppy 
Jo duj 



) 



Ap, (29) 

(30) 
(31) 



These generalized Manakov equations are deterministic, as 
they contain no rapidly fluctuating random terms. Equations 29 
represents an extension of the standard Manakov equations, 
for multimode fibers with birefringence. The first nonlinear 
term represents the intramodal nonlinear effects, resulting 
from self-phase modulation (SPM), and it occurs for single- 
mode fibers as well with the same coefficient 8/9. The second 
nonlinear term is new. Its origin lies in the intermodal nonlin- 
earities resulting from cross-phase modulation (XPM) among 
various fiber spatial modes. The averaging over birefringence 
fluctuations reduces the factor of 2 that is associated with XPM 
to a lower value of 4/3. Note that the random quantity q^p 
averages to zero. 

The new Manakov equations (29) obtained for multimode 
fibers can be solved numerically much faster than (12), and 
it should help considerably in designing future SDM systems. 
However, before it can be used with confidence, we need to 
ensure that its predictions agree with those of (12). In the next 
two sections we verify this by considering a specific SDM 
system. 

IV. Numerical simulations 

In this section we compare the performance of a specific 
SDM system for multimode fibers with and without random 
birefringence. Propagation equations are solved numerically 
with the split-step Fourier method. The case with no bire- 
fringence is obtained as a reference by solving (6). The case 
with birefringence is obtained either by solving explicitly (12) 
applying new random matrices at each step [20] or by solving 
the new set of generalized Manakov Equations (29). The 
comparison of the numerical results represents a test of the 
new Manakov equations. 

For simplicity, we assume that each spatial mode carries 
a single-wavelength channel (no WDM) at a bit rate of 
114 Gb/s. We make use of polarization multiplexing (PDM) 
together with the QPSK format so that the symbol rate is 
Rs = 28.5 Gbaud. Each PDM-QPSK signal is constructed 
from an independent random bit stream that is Gray-mapped 
onto the QPSK symbols. The pulse shape corresponds to a 
raised-cosine spectrum with a roll-off factor of 0.2. Each 
PDM-QPSK signal is made of 2^^ modulated symbols with 
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TABLE I: General shape of spatial mode distributions sup- 
ported by the step-index and graded-index fibers used in 
numerical simulations. 



Spatial distribution shapes 














Mm 





Step-index fiber 



LPOl 



LPlla, LPllb 



LP02 



LP21a, LP21b 



Graded-index fiber 



HGOO 



HGOl, HGIO 



HG02-hHG20 



HGlla, HGllb 



2^^ samples per polarization. The injected power per spatial 
mode is Pin = 7 dBm. This choice is deliberately high in order 
to observe large nonlinear impairments so as to represent a 
very stringent test of the new Manakov equations. . 

The 1000-km transmission line consists of 10 sections of 
100-km multimode fiber, each section being followed by an 
ideal erbium-doped fiber amplifier (EDFA) that can com- 
pensate for the span losses of 0.2 dB/km (assumed to be 
the same for each spatial mode). At the end of the fiber 
link, we artificially degrade the signal by adding noise to 
account for the amplified spontaneous emission (ASE) added 
by all EDFAs. This is justified since the nonlinear interaction 
between ASE and the signal has been verified to be much 
smaller than the nonlinear interaction of the signal on itself. 

We focus on a step-index multimode fiber with a core diam- 
eter of 12.3 /im and a numerical aperture of 0.2 (A = 0.01). 
Such a fiber has the V parameter of 5 at 1550 nm, and 
it supports LPOl, LP02, LPll, LP21 modes, resulting in a 
total of 6 spatial modes when we take into account 2-fold 
degeneracy of the LPl 1 and LP21 modes [17]. Table I presents 
the spatial distributions F^(x, y) of the different spatial modes 
supported by the fiber. 

Each spatial mode has its own propagation constant, modal 
group velocity and GVD parameter D (see Tab. II). The 
fundamental-mode nonlinear coefficient is 7 = 1.4 W~^km~^. 
We consider the case where there is no coupling between the 
spatial modes, i.e. qrnp = 0. 

We solve numerically (12) and (29) with a step size of 
100 m in all cases (larger step sizes could have been used 
in the Manakov case). For simplicity, we choose S$ip = in 
the case of (12). Smaller step sizes were tested and shown to 
give the same results. 

At the receiving end, the symbol stream is coherently 
detected by polarization diverse 90-degree hybrid followed 
by analog to digital converters. We assume an ideal DSP 
scheme that compensates perfectly all linear impairments 
such as chromatic dispersion, group-velocity mismatch, mode 



TABLE II: Differential modal group delay (compared to 
the fundamental mode) and group velocity dispersion D of 
each spatial mode of the step-index multimode fiber used in 
section IV 





DMGD [ns/km] 


D [ps/(km-nm)] 


LPOl 





25 


LPll 


6.5 


27.3 


LP02 


9.9 


-2.3 


LP21 


12 


20.8 



V LP01 

□ LP01 & LP11a 

O LP01 & LP11a& LP02 




16 18 
OSNR [dB] 

Fig. 1 : BER averaged over the co-propagating spatial modes of a step 
index multimode fiber for 1000 km transmission versus OSNR. 114- 
Gb/s PDM-QPSK signals are transmitted over different combinations 
of spatial modes indicated in the upper right inset. Solid curves 
obtained solving Manakov Eq. (29) are superposed with dashed 
curves representing the case with random birefringence obtained 
solving (12); pointed curves represent the zero-birefringence case (6). 
Back- to-back case (no fiber) is displayed as a reference. 

mixing, and polarization mixing. This is justified in practice 
since modern equalizers can remove most, if not all, linear 
impairments [21]. It also lets us present our results in a way 
that is independent of the type of equalizer employed. 

Figure 1 shows the bit-error rate (BER) averaged over all 
the propagating spatial modes as a function of the optical 
signal-to-noise ratio (OSNR) [1] after 1000 km multimode 
fiber when 114-Gb/s bit streams are transmitted on different 
combinations of spatial modes. Clearly, the predictions of our 
new Manakov equation (29) agree very well with those of 
the full stochastic equation (12). As discussed later, the use 
of Manakov equations reduces considerably the computation 
time. To study the impact of random birefringence on system 
performance, we also consider the case of an ideal multimode 
fiber with no birefringence by solving Eq. (6). It is evident 
that birefringence actually improves system performance in 
all cases by lowering the BER. It should be stressed that 
this performance improvement occurs in practice only if 
the random polarization mixing resulting from birefringence 
can be efficiently compensated at the receiver. The system 
performance is easily understood if we recall that the averaging 
over rapidly varying birefringence leads to (29) in which the 
nonlinear coefficient associated with the SPM is reduced from 
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12 14 16 18 20 22 



OSNR [dB] 

Fig. 2: BER of each spatial mode versus OSNR after 1000 km 
transmission on step-index fiber with rapidly varying birefringence. 
Solid curves represent the case where 6 spatial modes co-propagate; 
dashed curves represent the cases where spatial modes propagate 
individually. 

1 to 8/9. Note that for the chosen configurations of transmitted 
spatial modes, inter-modal XPM effects are averaged out 
because of the large differential group delay between spatial 
modes [22]. The reason of the large performance degradation 
observed in Fig. 1 in the case of 3 co-propagating spatial 
modes (LP01+LP11+LP02) is related to a small magnitude 
of the dispersion parameter D of the LP02 spatial mode for 
the chosen fiber design parameters. Figure 2 presents the 
received BER of each spatial mode. In one case, the 6 spatial 
modes are transmitted together. In the other case, each spatial 
mode is transmitted individually in the multimode fiber. The 
performance degradation between these two cases is due to 
the inter-modal XPM. We observe that LP02 and LPOl spatial 
mode suffer very little from inter-modal XPM because of their 
large differential group delay compared to the other spatial 
modes. On the other hand, performance of both LP 11 and 
LP21 modes are largely degraded by the inter-modal XPM 
because the degenerated spatial modes (for instance LPlla 
and LP lib) propagate at the same group velocity. Note that 
two degenerated spatial mode perform the same. The nonlinear 
inter-modal XPM coefficient is/ii^,ii^,iij,,iij,=0.3 between 
LPlla and LPllb and /21,, 21^,215,21, = 1 between LP21a 
andLP21b. 

In a step-index multimode fiber, there is a spread of 
group velocities between spatial modes. This group- velocity 
mismatch introduces linear degradations (similar to PMD) 
that must be compensated at the receiver. In practice, the 
equalization complexity grows in proportion to the magnitude 
of the differential group delays among modes. This problem 
can be solved using a graded-index fiber for which group- 
velocity mismatch nearly vanishes. It is thus important to study 
how the system performance in Fig. 1 will change if the step- 
index multimode fiber used there were replaced with a graded- 
index multimode fiber. In particular, we focus on a multimode 
fiber whose core has a quadratic graded refractive index (in 
the shape of a parabola). Such fibers are designed to have all 
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Fig. 3: BER averaged over the co-propagating spatial modes of a step 
index multimode fiber for 1000 km transmission versus OSNR. 114- 
Gb/s PDM-QPSK signals are transmitted over different combinations 
of spatial modes indicated in the upper right inset. Solid curves 
obtained solving Manakov Eq. (29) are superposed with dashed 
curves representing the case with random birefringence obtained 
solving (12); pointed curves represent the zero-birefringence case (6). 
Back- to-back case (no fiber) is displayed as a reference. 

spatial modes propagate at the same group velocity [17]. 

For numerical simulations, we consider a graded-index fiber 
with a core diameter of 17.4 /im, D = 21.5 ps/(km-nm) 
and 7 = 1.4 W~^/km and designed such that it supports the 
(Hermite-Gauss) modes presented in Table I. Data transmitted 
on different combinations of these spatial modes and we 
assume no linear coupling between them, i.e. q^p = 0. 
For similarity with LP02 spatial mode of step-index fibers, 
we consider a spatial distribution corresponding to the sum 
of the HG02 and HG20 spatial mode distributions. Figure 
3 shows the BER as a function of OSNR under the same 
conditions used in Fig. 1 except that the graded-index fiber 
replaces the step-index fiber. A comparison of two figures 
shows several interesting features. First, a rapidly varying 
birefringence results in a larger improvement in the case of 
graded-index fibers than in the case of step-index fibers. Here, 
spatial modes have the same group velocity so, there are 
nonlinear interactions between them due to inter-modal XPM. 
The effect of rapidly varying birefringence results to reduction 
of the XPM coefficient from 2 to 4/3, which improves the 
performance. Second, nonlinear penalties, for graded-index 
fibers, are larger in the two-mode case but become consid- 
erably smaller in the three-mode case. These results suggest 
that graded-index multimode fibers may perform better when 
a large number of spatial modes are used to simultaneously 
transmit data. 

Finally, we come to a main advantage of using the new 
Manakov equations (29). We have already seen that its nu- 
merical predictions are virtually identical to those obtained 
by solving Eq. (12). From a practical standpoint, the use of 
Manakov equations is preferable because it requires much 
less computational time. One reason is that the number of 
nonlinear terms is drastically reduced. For instance in the 
case of 3 propagating modes, there are 81 nonlinear terms in 



NONLINEAR PROPAGATION IN MULTIMODE FIBERS 



7 



TABLE III: Computing time in seconds for 1000-km transmission 
of 114-GB/s channels over multiple modes with a 100-m step size. 



Number of modes, M 


1 


2 


3 


4 


Manakov Eq. (29) 


270 


800 


1600 


2800 


unaveraged Eq. (12) 


720 


7600 


16000 


31000 



Eq. (12) but only 3 terms are needed in Eq. (29). In general, 
there are — 1 intermodal nonlinear terms for a M-mode 
fiber while there are M — 1 intermodal nonlinear terms for 
each generalized Manakov equation. A second reason is that 
a solution of Eq. (12) requires computation of an exponential 
matrix of size 2M x 2M, M being the number of spatial 
modes, which is a time-consuming numerical operation. As 
a result, the reduction in computing time depends on the 
value of M. This is evident in Table III, which shows the 
computation time obtained on a desktop computer using the 
MATLAB software (version 2007b) for 1000-km propagation 
of 114-Gb/s channels with a step size of 100 m and a temporal 
grid size of 2^^. Computation time is reduced by nearly a 
factor of 10 even for M = 2, with larger reductions occurring 
for larger values of M. Identical step sizes were used in 
both cases. Solving Eq. (12) numerically requires a step size 
that is determined by birefringence fluctuations. In the case 
of Manakov equations (29), the step size is mainly set by 
nonlinear effects which usually occur on a length scale much 
larger than that of birefringence fluctuations. Consequently it 
is usually possible to choose a much larger step size. Taking 
this into account, the use of the new Manakov equations should 
reduce the computing time easily by more than a factor of 100 
in most cases of practical interest. 

V. MuLTicoRE Fibers 

Manakov equations (29) in Section III assume weak cou- 
pling among distinct spatial modes of a multimode fiber. In 
practice, linear mode coupling can stem from such effects as 
core irregularities [23] or fiber bending [25]. Linear coupling 
between distinct spatial modes can be relatively weak in 
comparison to the coupling between the two polarization 
components of each mode. In our theory, linear mode coupling 
resulting from fiber imperfections is governed by the parameter 
Qmp{z) appearing in Eq. (6). This coupling coefficient varies 
randomly along the propagation distance and its strength 
depends on the degree of fiber imperfections and the level 
of stress applied to the fiber. 

Future SDM systems may also make use of multicore fibers 
in which each core supports a single mode, but modes in 
different cores can couple strongly if these cores are relatively 
close to each other. In this case, the extent of coupling depends 
on the physical distance among its cores and its magnitude 
varies exponentially with this distance. The coupling coeffi- 
cient qmp{z) may vary with z in practice due to waveguide 
imperfections but in order to simplify the analysis, we assume 
ideal fiber so that Qmp is a constant in the multicore study. 

For numerical simulations, we consider an ideal 3 -core fiber 
having identical, and equidistant single-mode cores, i.e., in 
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Fig. 4: BER averaged over the cores as a function of coupling 
length Lc for 1000 km transmission. 114-Gb/s PDM-QPSK with 
independent random bit streams are transmitted over a single-core 
(solid horizontal line) and a 3 -core (solid curve) fiber. Dashed curves 
show these two cases for fibers without birefringence. 

an equilateral triangular configuration. There are three spatial 
modes each associated with the field propagating in one 
core. The nonlinear propagation equation is still given by (6) 
but the coupling coefficients are now given by 

2Mml,Y'' II -"-)F™F,da:dy, (32) 

where nm{x^y) is the spatial distribution of the refractive 
index of the isolated cores in which the spatial mode m prop- 
agates. This propagation equation is obtained using a coupled- 
core approach, which assumes that the spatial distribution 
of a field propagating in one core, is not perturbed by the 
presence of other cores [23]. Because of the symmetry of the 
three-core fiber considered here, all coupling coefficients are 
identical and we can replace Qmp with a single parameter q. 
Each core of the fiber has a core diameter of 8.2 /im, with 
A = 0.003 and 7 = 1.4 W-^km-^ With these design 
parameters, intermodal nonlinear coefficients are negligible 
[24]. We modify the strength of linear coupling by adjusting 
the distance between fiber cores. It is useful to present the 
results in terms of the coupling length Lc = 7r/(3g), which 
represents the propagation distance at which maximum linear 
power transfer occurs among the cores. The rapidly varying 
birefringence is included by solving the propagation Eq. (12) 
with a step size of 100 m and changing the rotation matrix 
randomly after each step. 

Figure 4 shows the BER as a function of the coupling 
length (solid blue curve) when 114-Gb/s independent random 
bit streams are simultaneously transmitted through all three 
cores. The case of a single-core fiber is shown for comparison 
by a solid black line. Also shown for comparison are the two 
dashed curves for which the fiber link has no birefringence. 
The injected power per core is Pin = 8 dBm and the 
OSNR at the receiver is set to 20 dB in all cases. When 
the coupling length is comparable to the length scale Lf 
of birefringence fluctuations (here 0.1 km), there is a large 
performance improvement in the case of the 3 -core fiber. 
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This improvement can be understood by noting that a rapid 
power transfer among various modes induced by a strong 
Hnear coupHng and birefringence-induced polarization mixing 
helps to mitigate nonlinear effects. As predicted by the new 
Manakov equations (29) (where the term c^mp in (15) averages 
to zero), when the coupling length becomes much larger 
than Lf, the effects of linear coupling are averaged out by 
rapidly varying birefringence, and the 3-core fiber behaves 
similarly to a single-core fiber. This occurs in Figure 4 for 
Lc > 10 km. Indeed, when Eq. (29) is used, the numerical 
results correspond to the case of a single-core fiber. The main 
message from Figure 4 is that the new Manakov equations 
(29) can be used with success for multicore fibers as long 
as the cores are far enough apart that the coupling length is 
about 100 times larger than the length scale of birefringence 
fluctuations. 

These results also suggest that in the case of multimode 
fibers, Manakov equations (29) are valid when the coupling be- 
tween spatial modes, resulting from fiber imperfections, occurs 
on a length scale 100 times larger than that of birefringence 
fluctuations. 
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Fig. 5: BER averaged over the co-propagating modes versus OSNR 
after 1000 km over graded-index fibers. 114-Gb/s bit PDM-QPSK 
signals are transmitted on one, three and six spatial modes. Square 
and diamond markers represent the weakly and strong coupling 
regimes, respectively. The dotted curve shows the back-to-back case 
(no fiber). 



VL Manakov Equations in strong-Coupling 
Regime 

The Manakov equations (29) obtained in the weak-coupling 
regime fail if linear coupling among various spatial modes 
becomes comparable to birefringence-induced coupling. In 
this section, we derive new Manakov equations in the case 
of strong random coupling among M spatial modes of a 
multimode or multicore fiber. In this strong-coupling regime, 
we need to use 2M x2M random unitary matrices to account 
for random coupling among M spatial modes, with two 
polarization states per mode. 

Before applying such matrices to the nonlinear propagation 
equation, we need to rewrite the set (6) of M equations (p = 
1 . . . M) in a vector form. By introducing A = [Aj . . . Aj[^]^ 
as the column vector containing 2M field envelopes, we obtain 

dA dA B2 d'A 

^ [iATA)T*A* +2{A^J^A)J'A] dxdy 

(33) 

where J^isa2Mx2M matrix whose diagonal elements are 
M diagonal 2x2 matrices such that 



-Yi Fj I2, 



1 . . . M. 



(34) 



So, Bi and S2 are 2M x 2M diagonal matrices containing re- 
spectively the propagation constant, the inverse group velocity, 
and the dispersion parameter of each mode. Further, 



^So=So-^Tr(So) 



(35) 



represents deviations of the propagation constants from the 
average value of all spatial modes. Note that the coupling 
term qrap does not appear in (33) as it is included through 
the random unitary matrix IZ. 



Following the procedure in section III, we apply the substi- 
tution 

A = nA , (36) 

where 7^ is a random rotation matrix of dimension 2M x 2M 
whose elements rij satisfy the following relations [27]: 



/I |2| |2\ ^ ~^ ^kk 

{\rik\ \rik'\ ) 



{\rkmf\rk^m^f) 



{^ki'^kj^k'i'^k'j)'- 



2M(2M + 1) 

1 



(37) 



(2M + 1)(2M^ '''' 



-1 



(2M + 1)2M(2M- 1) 



k 7^ k'(39) 



By applying the substitution (36) and averaging the propa- 
gation equation (33), we obtain the following generalization 
of the Manakov equations in the strong-coupling regime: 



dA IdA [32 A ,,,2, 



where 



32 



Ikkll 



M 

^ 2^6M(2M + 1) ' 



(40) 



(41) 



1/v = trace (Si) /2M is the average inverse group velocity, 
and ^2 = trace (S2)/2M is the average group velocity 
dispersion. Equations (40) and (41) are in agreement with the 
results obtained in [16] using a different procedure. 

A comparison of the magnitudes of various nonlinear 
coefficients appearing in (29) and (40) reveals that system 
performance are improved in the high-coupling regime as the 
number of spatial modes increases. The nonlinear coefficient 
K, in Eq. (41) tends to zero with M increasing, which predicts 
that multimode fibers in the high coupling regime can perform 
better than single-mode fibers for large number of propagating 
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spatial modes. Using the same parameter values used in 
section IV, Figure 5 presents numerical results for 1000-km 
transmission over a graded-index fibers supporting one, three 
and seven Hermite-Gaussian HGmn spatial modes. Each fiber 
is designed in order to have 7 = 1.4 W~^/km. With the 
three mode fiber, data are transmitted on the HGOO, HGOl 
and HGIO spatial modes. With the seven mode fiber, data are 
transmitted on the six spatial modes presented in Table. I. 
Comparison of the weak and strong coupling regime shows 
that, as expected, the BER curves in the two cases coincide 
when only one spatial mode of the fiber is used. However, 
the high-coupling regime results in better BER performance 
when data is transmitted using multiple spatial modes. Another 
difference is that the performance is degraded in the weak- 
coupling regime as the number of spatial modes increases 
whereas, somewhat surprisingly, performance is nearly the 
same for three or six spatial modes in the high-coupling 
regime. 

VII. Conclusions 

In this paper we have focused on SDM systems designed by 
using multimode or multicore fibers and derived a nonlinear 
propagation equation satisfied by the bit stream transmitted 
through each optical mode in the presence of fiber disper- 
sion, birefringence, and nonlinearity. We first expressed this 
equation in Jones vector form in the slowly varying envelope 
approximation so that both polarization components of each 
mode can be treated simultaneously and we then used it to 
investigate the performance of SDM systems designed with 
multimode fibers exhibiting modal birefringence that varies 
randomly on a length scale smaller than the length over which 
nonlinear effects become important. The effects of fluctuating 
birefringence were included through random Jones matrices. 
As expected, the XPM and FWM effects resulting from 
intermodal nonlinearities cannot be ignored for multimode 
fibers. There are as many as — 1 intermodal nonlinear 
terms for an M-mode fiber. 

Since numerical simulations based on such stochastic non- 
linear propagation equations are generally time-consuming, we 
generalized the Manakov equations, well known in the case of 
single-mode fibers, to the case of multimode fibers. We con- 
sidered the weak linear coupling case first and averaged over 
birefringence fluctuations following a procedure similar to that 
used in the original derivation of the Manakov equations for 
single-mode fibers. We found that averaging over birefringence 
fluctuations reduces the — 1 intermodal nonlinear terms to 
only M — 1 terms because all FWM-type terms disappear. 
The XPM-type terms remain but their effectiveness is reduced 
because the well-known factor of 2 is reduced to 4/3. 

We verified the validity of our new Manakov equations 
by simultaneously transmitting multiple 114-Gb/s bit streams 
in the PDM-QPSK format over different modes of a multi- 
mode fiber and comparing the numerical results with those 
obtained by solving the full stochastic equation explicitly. 
The agreement between the two methods was excellent in all 
cases studied. We found that the use of the new Manakov 
equations can reduce computation time by at least a factor 



of 10 and by a factor in excess of 100 as the number of 
modes grows. Our numerical results show that birefringence 
fluctuations improve system performance by reducing the 
impact of fiber nonlinearities. The extent of improvement 
depends on the fiber design and how many modes are used 
for SDM transmission. We investigated both the step-index 
and graded-index multimode fibers and considered up to three 
spatial modes in each case. 

We extended our theory to the case of multicore fibers 
and discussed the validity of the new Manakov equations by 
studying the impact of linear coupling among various cores. 
We found that the new Manakov equations remain valid when 
the length scale of random birefringence of the cores is about 
100 times smaller than the linear coupling length among the 
cores of a multicore fiber. 

Finally we studied the case of strong random coupling 
among various spatial modes of a multimode fiber and derive 
new Manakov equations for this regime. We showed though 
numerical simulations that the high-coupling regime leads to 
better system performance than the weak-coupling regime 
when the number of spatial modes is at least two. Theory 
also predicts that performance improves for large number of 
propagating spatial modes so that such multimode fibers can 
perform even better than single-mode fiber. 
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